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■ Abstract. Newtonian point mass binaries can be brought into arbitrarily close circular 
\ orbits. Neutron stars and black holes, however, are extended, relativistic objects. Both 

finite size and relativistic effects make very close orbits unstable, so that there exists an 
innermost stable circular orbit (ISCO). We illustrate the physics of the ISCO in a simple 
model problem, and review different techniques which have been employed to locate 
the ISCO in black hole and neutron star binaries. We discuss different assumptions 
^ ■ and approximations, and speculate on how differences in the results may be explained 

lO ' and resolved. 

o 

O : INTRODUCTION 

o 

■ Compact binaries, containing black holes or neutron stars, are among the most 
O^' promising sources for the new generation of gravitational wave detectors. TAMA 

! has already started taking data, and LIGO, VIRGO and GEO will soon become 

^1 operational (see, e.g., [1]). With the advent of gravitational wave astronomy arises 

a need for theoretical templates of gravitational waveforms, which are required for 
^ ■ the identification and interpretation of signals in the noisy output of the detectors. 

Compact binaries emit gravitational radiation, and therefore loose energy and 
slowly spiral towards each other. Because of the circularizing effects of gravitational 
radiation damping, it is reasonable to assume the orbits of close binaries to be 
quasi-circular. The slow and adiabatic inspiral continues until the binary reaches 
the innermost stable circular orbit (ISCO), at which the orbits become unstable. 
At that point, the stars start to plunge towards each other, and coalesce and merge 
after a dynamical timescale. The ISCO leaves a characteristic signature in the 
gravitational wave signal of a binary inspiral, and therefore provides an important 
piece of information for the construction of gravitational wave templates. 

In this article, we will explain in a simple point-mass model problem how both 
tidal and relativistic effects make very close orbits unstable, giving rise to the ISCO. 
We will then review attempts to determine the ISCO for binary black holes, which, 
for stellar-mass black holes, is likely to fall into the frequency range of LIGO. In 
particular, we will compare results from post-Newtonian (PN) and numerical calcu- 



lations. We will also discuss some of the qualitative results for binary neutron stars, 
and will summarize some of the more recent results from relativistic simulations. 



THE ISCO IN A SIMPLE MODEL PROBLEM 

Most commonly, the ISCO is located with the help of turning-point methods. We 
motivate this method by applying it to Newtonian point-masses, and introduce tidal 
and relativistic correction terms to illustrate their effect. Rigorous justifications for 
turning-point methods have been developed, for example, in [2-4]. 



Newtonian point- masses 



The conditions for a circular orbit can be derived very easily from the Hamilto- 
nian formalism. Consider, for example, a Newtonian point mass binary, for which 
the energy is the sum of the kinetic and potential energy 
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Here M = mi + m2 is the total mass, ji = niim2/M the reduced mass, r the 
separation vector, and — the angular velocity. Rewriting the energy in terms 
of the conjugate momenta P — /ir and J = /ir^fl yields the Hamiltonian 
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The four independent variables r, P, tp and J now satisfy the Hamiltonian equations 
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The last equation shows that the angular momentum is conserved, since the Hamil- 
tonian is cyclic in <p. To construct a circular orbit, we obviously need r = 0, which 
imphes P — 0, and also P — 0, so that 
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(4) 



A circular orbit hence corresponds to an extremum of the energy at constant angular 
momentum. The orbital frequency of a circular orbit can be determined from 
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FIGURE 1. Energy as function of separation for a Newtonian point mass binary (left panel), 
and a model problem including tidal effects (right panel, see text). The solid lines are contours of 
constant angular momentum J, with the values of J increasing from the bottom to the top. The 
extrema of these contours correspond to circular orbits (eq. (4)). The dashed line connects the 
circular orbits and represents the equilibrium energy Eeq, and the dots mark the ISCO for tidal, 
relativistic and combined effects. 

The last two equations are crucial for the construction of the ISCO. 

Returning to the example of a Newtonian point mass binary, we can construct 
circular orbits by inserting the energy (2) with P — into eqs. (4) and (5). Eq. (4) 
yields the virial relation 



so that the equilibrium energy of a circular orbit becomes 



The orbital frequency of these orbits can be found from eq. (5) and satisfies, as 
expected, Kepler's third law 
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where we have used eq. (6) in the last equality. 

Note that Newtonian point masses can be brought into arbitrarily close and 
arbitrarily tightly bound circular orbits. The orbits are stable for all separations r. 



indicating that in Newtonian point mass binaries there is no ISCO. As we will see, 
this is related to the gravitational interaction potential being proportional to 1/r. 

It is also illustrative to construct circular orbits graphically. In the left panel of 
Fig. 1, the solid lines are contours of the energy (2) at constant angular momentum. 
The minima of these contours correspond to circular orbits. Connecting these yields 
the equilibrium energy (7) as a function of separation. 



Tidal and relativistic effects 



Compact binaries are extended, relativistic objects, of course, so that we have to 
take into account both tidal and relativistic effects. 

For irrotational, identical stars, tidal effects can be modeled by adding a tidal 
interaction term 

£^tidai = -2A^^^ (9) 

to the energy (2) (compare [5-7]), where R is the radius of the stars at large 
separations. The coefficient A depends on how easily the stars can be deformed, 
and hence on their structure and equation of state (EOS). For polytropic EOSs 
with polytropic index n we have A = (3/4) «„(! — n/5), where the coefficient «;„ 
is tabulated, for example, in [2]. For an incompressible fluid, A = 3/4. Each star 
induces a quadrupole moment in the companion^ which scales with 1/r^, and the 
two quadrupole moments' interaction again scales with 1/r^, so that £^tidai ^ 
From eq. (4), we can now find the equilibrium energy 

We immediately see that we can no longer construct arbitrarily tightly bound orbits, 
and that instead E^q now assumes a minimum for a finite separation r. It is easy 
to show that any attractive interaction potential proportional to 1/r", n > 2, leads 
to a positive contribution to E^q, and hence the existence of a minimum at finite r. 

In the right panel of Fig. 1 we provide a graphical representation for incompress- 
ible fiuid stars with m/R = 0.2, where m is the mass of the individual stars. For 
large separations, the tidal interaction is very small, and the orbits are similar to 
those of point masses. For small separations, however, the tidal interaction becomes 
dominant, and decreases the energy of a configuration with a given separation and 
angular momentum below that of the point mass binary. For large enough angular 
momentum, the J = const contours now have a maximum at a small separation in 
addition to the minimum at a larger separation, while for small angular momenta 
their contours no longer have any extrema. As a consequence, the equilibrium en- 
ergy goes through a minimum. Outside of this minimum, the equilibrium energy 



^) If the stars are not irrotational, the spin of the individual stars induces an intrinsic quadrupole 
moment. The tidal interaction energy then scales with a smaller power of r, see, e.g. [5,25] 



curve connects minima of the J = const curves, which correspond to stable circular 
orbits, while inside this minimum, the curve connects maxima of J — const curves, 
which correspond to unstable circular orbits. The minimum of the equilibrium 
curve therefore marks the innermost stable circular orbit, the ISCO^. 

We can similarly mimic relativistic effects by borrowing a relativistic interaction 
potential 

Erel = T 11 

from the result for test particles in Schwarzschild space times. This interaction 
again gives rise to an ISCO. We mark the location of this ISCO in the left panel of 
Fig. 1, as well as the ISCO when computed from the combined tidal and relativistic 
effects. It is evident from this model calculation that including relativistic effects 
will move the ISCO to larger separation, and correspondingly smaller values of the 
orbital frequency 

Obviously, this model problem is very naive, and can at best mimic qualitative 
effects. However, it does illustrate several important results. Quite in general, 
the ISCO arises because of contributions to the interaction potential which deviate 
from a simple 1/r scaling. In particular, tidal effects can cause an ISCO even in 
purely Newtonian systems. In general, both tidal and relativistic effects have to be 
taken into account for an accurate determination of the ISCO. Lastly, the model 
problem provides a straight-forward recipe for locating the ISCO: first construct 
circular orbits and the corresponding equilibrium energy E^q by locating extrema 
of the energy (eq. (4)), then identify the ISCO with the minimum of E'eq; and 
last compute the orbital frequency fl from eq. (5). We have therefore reduced the 
problem of accurately determining the ISCO in black hole or neutron star binaries 
to an accurate determination of their energies, which we will address separately in 
the following sections. 



THE ISCO IN BINARY BLACK HOLES 

Determining the ISCO in binary black hole systems has been attempted with two 
independent approaches: post-Newtonian expansions and numerical calculations in 

At this point, a word of warning is in order. Relativistic binaries emit gravitational radiation, 
causing them to slowly spiral towards each other, and they hence do not follow strictly circular 
orbits. The very concept of an innermost stable circular orbit is therefore somewhat ill defined. 
Also, the minimum in the equilibrium energy identifies the onset of a secular instability, while 
the onset of dynamical instability may be more relevant for the binary inspiral (see, e.g., the 
discussion in [2] and also [3], where it is shown that the two instabilities coincide in irrotational 
binaries). Moreover, it has been suggested that the passage through the ISCO may proceed quite 
gradually [8], so that a precise definition of the ISCO may be less meaningful than the above 
turning method suggests. Ultimately, dynamical evolution calculations will have to simulate the 
approach to the ISCO and to investigate these issues. For the sake of dealing with a well-defined 
problem, we will here identify the ISCO with the minimum of the equilibrium energy. 



full general relativity. We will outline both approaches, and will then compare the 
results. 



Post-Newtonian Calculations 

A large number of researchers has developed post-Newtonian techniques to model 
compact objects in close binaries (an incomplete hst includes [9-13]; see [14] for a 
particularly pedagogical treatment) . Typically, these calculations start by bringing 
Einstein's equations into the form 

nh""^ = -lenr''^ = -167r(-^)T"^ - A"^. (12) 

Here, □ is the fiat space wave operator, h°'^ measures deviations from the fiat 
Minkowski metric, g is the determinant of the metric, and the source term r"^ 
contains both the stress-energy tensor T"^ as well as non-linear terms in h"'^, which 
have been absorbed in A"^. A formal solution to this equations is the retarded, 
flat-space Green function 



Unfortunately, the source term r°'^ depends on the solution h"'^, so that this formal 
solution is of little practical help. However, it can be used to construct a solution 
iteratively. Starting with a Newtonian point-mass solution, for which is 
known and for which h"^ vanishes, we can construct a first iteration from 

nhf = 16nTSL- (14) 

Each solution h^'^ can be inserted on the right hand side of eq. (13), which then 
provides the next iteration h^+i- Each iteration provides a correction over the 
previous one in the order of e ~ ~ m/i?. The iteration therefore yields a post- 
Newtonian expansion of the solution /i"^, from which the energy can be constructed 
in the form of a Taylor expansion 

^eq = ^Newt + 6 + 6^ + . . . . (15) 

Here E^ e" is the n-th order post-Newtonian correction to the energy. 

It turns out, however, that for values close to the ISCO, e ~ 1/6, this expansion 
converges extremely slowly. The reason is that at a location fairly close to the 
ISCO, namely at the light radius with e ~ 1/3, some of the functions involved in 
the expansion diverge and have a pole (see [10]). It is therefore to be expected 
that a polynomial expansion of the form (15) will converge only very poorly in 
the neighborhood of that pole. Damour, Iyer and Sathyaprakash [10] provided a 
solution to this problem by using a resummation technique. Instead of employing 
the polynomial (15), they use the information contained in the Taylor expansion 



to construct an expansion in terms of rational functions. To illustrate this with an 
example, assume that we know the Taylor expansion of an arbitrary function f{x) 
up to second order. We can then construct an expansion in terms of a rational 
function 

f{x)r^to + tix + t2X+...r^— ■ (16) 

l + qix + ... 

by choosing the coefficients po, pi and qi such that the value and the first two 
derivatives of the two expansions match at a; = 0. This expansion is known as a 
"Pade-approximant" , and has been shown to greatly improve the convergence of 
post-Newtonian expansions at least up to second order. 

When h^^ is inserted on the right hand side of eq. (13), the integrals no longer 
have compact support, so that there is no guarantee that they will converge. More- 
over, the use of point-mass sources gives rise to divergent integrals, which have to be 
re-normalized appropriately. This has been achieved up to second post-Newtonian 
order, but at third post- Newtonian order some ambiguities remain [11]. Damour, 
Jaranowski and Schafer [13] have recently shown that these ambiguities can be 
expressed by a single dimcnsionlcss parameter, ecstatic j which is currently unknown. 
They compare three difi:erent post-Newtonian approaches (the e-method and j- 
method based on minimizations of the energy and the angular momentum, and 
an effective one-body method [12]) and find a quite remarkable result. For one 
particular choice of cc^static, these three independent approaches yield very similar 
predictions for the angular momentum, energy, and orbital frequency at the ISCO. 
This self-consistency suggests that the correct values of a;static and the ISCO at 3PN 
order may have been identified. 



Numerical Calculations in General Relativity 

A framework for numerically constructing models of binary black holes has been 
provided by Arnowitt, Deser and Misner's 3-1-1 (ADM) decomposition of Einstein's 
equations [15] and York's conformal decomposition [16]. 

The ADM formahsm sphts Einstein's equations into constraint and evolution 
equations. The gravitational fields, described by the spatial metric '-/ij and the 
extrinsic curvature Kij, have to satisfy the constraint equations on each time slice, 
while the evolution equations describe how they evolve from one time-slice to the 
next. For the construction of initial data, the two constraint equations, the Hamil- 
tonian constraint and the momentum constraint, have to be solved, which determine 
only the longitudinal parts of the gravitational fields. The transverse parts of the 
fields, loosely associated with the gravitational wave degrees of freedom, are uncon- 
strained by the constraint equations, and have to be chosen before the constraints 
can be solved. 

Since the binary inspiral outside the ISCO proceeds very slowly, the gravitational 
wave content of these spacetimes must be very small. It therefore seems reasonable 



to attempt to minimize the gravitational wave content by choosing the spatial 
metric conformally flat, 7^- = i^^fij, where ^0 the conformal factor and fij a flat 
metric. With the further assumption of maximal slicing, K = — 0? ^^e 

Hamiltonian constraint reduces to 

= -l^|^-'AJA'^, (17) 

while the momentum constraint becomes 

ViA'^ = 0. (18) 

Here V is the flat space covariant derivative, and A^^ is the trace-free part of the 
conformally related extrinsic curvature, Aij = ijj'^{Kij — '~fijK/3). Quite remark- 
ably, the momentum constraint (18) is now a linear equation and decouples from 
the Hamiltonian constraint. Solutions describing a pair of black holes with arbi- 
trary momenta and spins can now be constructed analytically by super-imposing 
two solutions for single black holes [17]. Given these solutions, the Hamiltonian 
constraint (17) can then be solved numerically [18,19], which completes the con- 
struction of binary black hole initial data. 

To construct binary black hole models in circular orbits and to determine their 
ISCO, turning points of the energy of these solutions have to be located, in complete 
analogy to the model problem above [20.21]'^. These results have recently been 
generalized to binaries in which the individual black holes carry spin [22]. 



Comparison 

In Fig. 2, we summarize results from the post-Newtonian and numerical calcu- 
lations. It is obvious that both the 2PN and 3PN results differ quite significantly 
from the numerical results, by more than a factor of two in the frequency. This 
discrepancy is very unsatisfactory, and should be resolved by re-examining both 
approaches and their assumptions. 

The PN approach treats the binary stars as point-masses, and hence neglects tidal 
effects (even though these could be included, see, e.g.. Appendix F of [14]). This 
approximation is probably fairly poor for binary neutron stars, but may be more 
adequate for binary black holes. Moreover, including tidal effects would probably 
move the ISCO to a larger separation and smaller orbital frequency, and would 
hence increase the discrepancy between the PN and numerical results. However, 
the convergence properties of the PN expansion and the lacking renormalization of 
higher-order expansion coefficients remain worrisome. 

The two approaches of [20] and [21] differ in the topology chosen for the binary black hole: [20] 
adopts the two-sheeted topology of [18], whereas [21] assumes the three-sheeted topology of [19]. 
Their results agree to within a few percent, indicating that this choice has little effect on the 
ISCO. 
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FIGURE 2. Results for the angular momentum, energy and orbital frequency of a black hole 
binary at the ISCO from post-Nc-wtonian and numerical calculations. The solid triangle, square 
and circle are the 2PN results of [13] using the effective one-body method, the e-method, and the 
J -method. The open circle is their 3PN result, -with the unknown parameter Wstatic chosen such 
that all three methods agree. The four-pointed and six-pointed star arc the numerical results 
of [20] and [21]. The dashed triangle and square are 2PN results using the effective one-body and 
e-method together -with a conformal-flatness assumption. The top label gives the corresponding 
gravitational wave frequencies for a binary of two 25 Mq black holes. 



The biggest worry in the numerical calculations is probably the assumption of 
conformal flatness. The effect of this assumption is evaluated in Appendix B of 
Ref. [13], where the authors determine the ISCO at 2PN order, using their effec- 
tive one-body and e-method together with a conformal-flatness assumption. These 
results are included as the dashed triangle and square in Fig. 2. Obviously, this 
assumption changes the results somewhat, but not nearly enough to explain the 
discrepancy between the PN and numerical results. When evaluated for the effec- 
tive one-body method, it even seems like conformal flatness is quite an adequate 
approximation. 

None of the above more obvious worries seem to be able to completely explain 
the differences between the PN and numerical results. This suggests that perhaps 

the two approaches construct sequences which are not physically identical. In our 
discussion of the model problem above, we have implicitly assumed that the masses 



mi and m2 of the binary stars remain constant during the inspiral. In general 
relativity, however, there is no unique definition of the mass of the individual black 

holes in a binary, and it is not clear which mass is conserved during the inspiral. 
It has been conjectured [23] that during the adiabatic inspiral outside of the ISCO 
the "irreducible mass" [24] of the black hole event horizons is conserved. In the 
numerical calculations of [18,21], this is approximated by the irreducible mass of the 
black hole apparent horizons'^. In the PN approaches the masses mi and m2 of the 
point sources are kept constant. It is not obvious that the two approximations are 
equivalent, and it is therefore possible that the two approaches construct physically 
distinct sequences. 



THE ISCO IN BINARY NEUTRON STARS 

The ISCO in binary neutron stars has been computed in numerous Newto- 
nian [2,5,25-27], PN [3,7,28] and relativistic calculations [29-31]. We wiU briefly 
summarize some of the qualitative flndings of these calculations, and will then 
discuss some of the more recent relativistic results. 



Qualitative discussion 

For binary neutron stars, the location of the ISCO depends on the physical 
properties of the individual stars, in particular their EOS and spin. Qualitatively, 
these effects can be illustrated by re-examining the tidal correction to the energy (9). 
Computing Eeq from (4) and (9), and then locating the minimum of Eeq, shows 
that the ISCO occurs at a separation 

^ = (48A)V5 (19) 
R 

and a frequency 

mQisco^yf (48A)-^A°g)'''. (20) 

Note that risco scales with the stellar radius R, and mf^isco accordingly with 
(m/i?)^/^, where m/R is the stellar compaction (compare [25,7,31]). 

The coefficient A depends on the EOS. For polytropic EOSs, A increases with 
decreasing polytropic index n. Eq. (19) therefore implies that risco/R is larger 
for stiffer EOSs. This explains why an ISCO exists only for stars with sufficiently 
stiff EOSs; for stars with too soft EOSs the stars merge before they encounter the 
ISCO. The critical value of n for the existence of an ISCO depends on the rotation 



Locating an event horizon requires knowledge of a complete spacetime, while an apparent 
horizon can be located on a single timeslice. 



of the stars and on relativistic effects, but is generally believed to be fairly close to 
JT-crit ~ 1- Eq. (20) also implies that the ISCO frequency is correspondingly smaller 
for stiffer EOSs. 

In the above discussion we have assumed that the stars are irrotational, so that 
Etidai ~ r^^. The effect of individual rotation can be estimated by allowing E'udai 
to scale with a different power of r (compare [3,5,25]). 

Numerical Calculations in General Relativity 

Constructing relativistic binary neutron stars requires solving the equations of 
relativistic hydrodynamics together with the initial value equations of general rela- 
tivity. This problem can be simplified whenever the fiuid fiow is stationary, so that 
the equations of hydrodynamics reduce to a relativistic Bernoulli equation^. 

This approach was first adopted by [29], who constructed corotational models of 
binary neutron stars in quasi-equilibrium. The initial value equations of general 
relativity were solved with the assumption that the spatial metric is conformally 
fiat. For moderate (but reahstic) compactions {m/R ^ 0.2), for which the gravi- 
tational fields are only moderately strong, the latter assumption has been shown 
to be quite adequate (compare [32]). Results for various polytropic indices n and 
stellar compactions m/R can be found in [29]. For stars of 1.4 Mq, n = 1 and 
m/R — 0.2, for example, they find a gravitational wave frequency at the ISCO of 
Qpsco = 2Qisco~ 1300Hz. 

The viscosity in neutron star interiors is not believed to be strong enough to 
maintain corotation during binary inspiral [33]. It is therefore more realistic to 
assume the binary to be irrotational, in which case the relativistic equations of 
hydrodynamics can again be reduced to a Bernoulli equation [34]. Constructing 
irrotational binaries is computationally more complicated than constructing coro- 
tational binaries, because one of the boundary conditions has to be imposed on 
the surface of the stars. The location of the latter changes during the numerical 
iteration, and is a priori unknown. Nevertheless, several groups have succeeded in 
constructing relativistic, irrotational models of binary neutron stars [30] . 

Probably the most careful analysis to date of evolutionary sequences of irrota- 
tional binary neutron stars is presented in [31]. As pointed out earlier [30], the 
authors find that for moderately soft EOSs (n > 2/3), a cusp forms at the surface 
of the stars before an ISCO is encountered, at which point the numerical method 
breaks down. Most likely, the cusp indicates that a Lagrange point forms, and that 
matter starts to overfiow towards the companion. The assumption of a stationary, 
irrotational fiuid fiow seems no longer adequate, and will have to be relaxed for 
the construction of closer binaries. An ISCO appears in these simulations only for 
fairly stiff EOSs {n ^ 2/3). For these cases, the authors of [31] quite carefully 

Note that a solution to the BernoulU equation is by construction in equiUbrium, and yields 
the equilibrium energy E^q. For binary black holes the latter has to be constructed by finding 
turning points of the energy E. 



analyze the ISCO for various compactions m/R, and show that it is dominated by 
the hydro dynamical effects of the tidal interaction. Up to moderate compactions 
{m/R < 0.2), the relativistic effects can be expressed as PN corrections to the 
Newtonian scaling (20). 

SUMMARY 

Knowledge of the ISCO may be an important ingredient for the future detection 
of gravitational wave signals and their interpretation. In this article we illustrate 
in a simple model problem how turning-point methods can be used to locate the 
ISCO (see, however, the disclaimer in footnote (2)). We also summarize recent 
efforts to determine the ISCO in both black hole and neutron star binary systems. 

Currently, the biggest challenge for an accurate determination of the ISCO in 
binary black hole systems seems to be the disagreement between PN and numer- 
ical approaches. For binary neutron stars, wc seem to be lacking a realistic de- 
scription of the velocity field in close binaries. Ultimately, the quasi-equilibrium 
approaches presented in this paper should be backed up with fully dynamical sim- 
ulations (e.g. [35]). 

As a final comment, we point out that in this article we have only discussed bi- 
naries containing either black holes or neutron stars. So far black hole-neutron star 
binaries have only been modeled in Newtonian or PN dynamical simulations [36] 
and in ellipsoidal model calculations [37]. Such mixed binaries have yet to be 
constructed in a relativistic framework, and their ISCO (or the onset of tidal dis- 
ruption) has yet to be determined self-consistently. 
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